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An attemr-t has been made to improve the gel and 

glass effect correlations of Chiu et al. Soong and 

coworkers adopted the model of Chiu et al. and used 

different values of the parameters in the correlations 

f'or and 9 . Ko consistent set of ^■'aremeters v^ere 
t F' 

obtained in their studies. Improved parameter values 
obatined herein match the CSTR simulation results of 
Carrat et al. very well. Steady state multiplicity 
patterns of continuous flow stirred tank reactors (CSTRs) 
V7ith gel effect are explored. Above an average conversion 
of 40 %, the p'resence of gel effect cannot be neglected. 
Only hysteresis type multipdicity could be obtained using 
simulation results. The stiffness and sensitivity of the 
defining equations of hysteresis and isola varieties 
made it impossible to solve them numericslly to give a 
bifurcation diagram. 



CHAPTER I 


INTRODUCTION 

Itiere has been a spate of research activity in the 
area of modeling and control of free radical polymerization 
reactors over the past few years. With the improved under- 
standing of polymerization kinetics and reactor behaviour, 
aided by the implementation of sophisticated reactor control 
algorithms, current research is directed towards the control 
of polymer properties. For over four decades, both free 
radical bulk and solution polymerizations have been a matter 
of study. From the vast array of experimental data with 
differoit systems (e.g. , polystyrene, PS, polymethyl methacrylate,! 

PMI'4A) in this field, it is well established that the classical 

1 2 

low conversion free radical kinetics ' (See Table I for 

mechanism) cannot be applied over the entire course of the I 

reaction, especially at high conversions. Even under I 

isothermal conditions, many such polymerization systems 

asdiibit a large auto acceleration of the rate of polymerization, 

associated with an increase in the molecular wei i^t of the 

polymer produced. Ihis auto acceleration, the result of a 

very rapid decrease in the rate constant associated with chain- 

termination (k^ and k^,) due to the introduction of severe 
tc td 

diffusional limitations, is known as the gel or Trommsdorf f 
effect. With the progress of reaction beyond the gel-effect. 
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TABLE I 


KINETICS OF FREE RADICAL POLYMERIZATION 


INITIATION 


k 


2R- 


PROPAGATION 


R. + M >.P. 


P .. + M 
3 




TERMINATION 


^td 


P. . + P .. 
1 3 




by disproportionation 




"tc 


I D. , . by combination 

^ i + j 
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a point may be reached where even the diffusion of the 
monomer to the reactive radical sites is severely curtailed. 
This phenomenon accounts for a drastic decrease in the 
propagation rate constant, k , and is known as the glass 
effect. In the past few decades , several researchers have 
tried to explain experimental data by improving the kinetic 
model to account for the gel and glass-effects. All such 
improvements are backed by some physical assijimptions. Friis 

4 5 

and Hamielec and Ftoss and Laurence used simple empirical 

methods to describe such polymerization kinetics. They 

correlated an apparent termination rate constant with various 

system parameters, such as conversion, temperature and free 

volume. A second string of workers approached the p3X»blem 

of gel effect by making some empirical statements about the 

molecular wei^t and conversion dependence of the termination 

rate constant, k^, using either qualitative analogies or 

experimental results on polymer transport phenomena or both. 

Of these workers, the work of Cardenas and O 'Driscoll^ (COD) 

.7 

and Marten and Hamielec (MH) are important. In the COD 

model, the polymer species are divided into two populations, 

one with the degree of polymerization (DP) less than a 

critical value and the other with a degree of polymerization 

greater than the critical value. The critical degree of 

polymerization dictates the onset of the gel effect. The 

population with a DP less than (DP) is assigned a 

termination rate constant 'k. ' which is the value in the 

: r> ■ 
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absence of the gel-effect. But the termination rate 
constant for the population with DP > DP^rit^ ^t' 
apparent rate constant, is an inverse function of the 
number-average degree of polymerization and the polymer 
volume fraction. In the MH model it is assumed that 
depends on the polymer concentration and molecular weight. 

Both these model are good in fitting experimental data, 
but are bsed on totally different and nonequivalent physical 
assumptions and use different dependences of upon 
concentration and molecular weight. However, the above 
models are empirical and are not tested against other 
solution properties of the syston over the range of 
conditionsof the experiment. Moreover, in these models, 
the diffusion behaviour of the radicals has been expressed 
in an adhoc manner which is flexible enough so that it 

could be adjusted to fit experimental data. Tulig and 

8 

Tirrel developed a model using classical diffusion- 
controlled reaction kinetics and related the termination 
rate constant to the diffusion constant of the growing 
radical. "Ihe model can fit experimental data on conversion 
vs. time and molecular weight vs. conversion upto about 

9-12 

70 or 75% conversion. Besides these, a number of researchers 
in an attempt to determine the effects of concentration, 
viscosity and molecular wei^t on the rate of polymerization, 
have carried out kinetic experiments. But even with such 
extensive research, a consistent understanding of the 
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molecular processes controlling the rate of polymerization 

13 

could not emerge. Chiu et al . recently focussed on this 
aspect and modeled this phenomena using a molecular view 
point. Ihey developed a mathematical model to describe the 
gel and glass effects from a fundamental molecular phenomena. 

It does not need any artificially imposed breaJc points either 
in the viscosity or in the molecular weight vs. concentration 
plotS/ to indicate the onset of diffusional limitations to 
the termination and propagation steps. Model parameters 
are linked to the operating variables of the process^ especially 
to the reaction temperature. Moreover, one unified equation 
is sufficieit to evaluate k and k. for the aitire course of 

p t 

the reaction. They have fitted the experimental data at all 

the three temperatures, 50°, 70° and 90°C, for which good 

isothermal batch data is available, using gross values of the 

model parameters. For these reasons , the model of Chiu 
13 

et al. is followed in the present work. 

In the last two decades, ranarkable progress has also 
been made in the understanding of the complex steady-state 
multiplicity patterns shown by continuous flow stirred tank 
reactors (CSTRs). For a single, exothermic, first order 
reaction occuring in a CSTR, changes in the flow rate (hence, 
residence time) can lead to complex multiplicity of steady- 

14 ^ ' 

states of the reactor. Zeldovich and Zysin were the first 
to indicate four diff erent multiplicity 'patterns for a very 
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high activation energy of the reaction. Although earlier 

unnoticed, this work was later confirmed by Furusawa and 

15 I 5 17 

Nishimura , and Hlavacek et al. . Uppal , Ray and Poore 

developed an intricate scheme to investigate the multiplicity 

pattern for any specific set of parameters, which was later 

18 19 

simplified by Huang and Varma . Golubitsky and Keyfitz 

applied singularity theory to prove that, for such systems, 

upto six different steady-state multiplicity patterns are 

possible (Figure 1). However, their analysis cannot predict 

the multiplicity pattern for any set of physical parameters 

Vcilues, Balakotaiah and Luss^*^ developed a simple technique 

for the exact prediction of the multiplicity pattern for any 

set a parameter values, ihey extended their technique to 

find the multiplicity pattern for reactions with any rate 
21—25 

expression . Althou^ these are devoted only to single, 

exothermic reactions, considerable progress has been made 
in its application to polymerization reactors also, indeed, 
in recent years, steady state characteristics and dynamic 
behaviour of various types of polymerization reactors have 
been studied extensively , Hamer et al. studied the 

free radical solution polymerization in CSTRs and identified 
multiple steady states, limit cycles and oscillatory pheno- 
mena depending upon the values of the solvent volume fraGtion 
and the heat transfer parameter. They showed that upto six 
differe:it types of multiplicity patterns are possible. 
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34 

Schmidt et al. confirmed the isolated branch in the 

conversion vs, residence time plot (called isola) 

experimentally. They indicated that operating the reactor 

at isolas may help obtain hi^ conversions at low residence 

times, ihey used coolant temperature, reactor heat transfer 

coefficient and feed solvent concentration as parameters 

and showed that these parameters can have an enormous affect 

28 

on reactor behaviour, jaisinghani and Ray performed an 

analysis of the reactor bdaaviour by assuming that the 

initiator concentration remain unchanged for the course 

of the reaction i.e., initiator consumption is negligibly 

small. They argued that, in commercial practice, the 

effective initiator concentration in the reactor can be 

maintained nearly constant over a wide range of conversions, 

when the initiators are slowly decomposing and are used in 

hi^ concentration in the feed. They identified the multiple 

35 

steady states and unstable limit cycles. Recently, Choi 
followed the technique of Balakotaiah and Luss^*^ to investi- 
gate the mtil tiplicity patterns of CSTRs for free radical 
solution polymerization. He developed the model for solution 

polymerization (solvent volume fraction f _ > 0,6) of methyl- 

s 

methacrylate with constant initiator concentration and no 
gel effect and analysed it using heat of reaction, heat 
transfer coefficient, coolant temperature and activation 
energies in the form of various dimensionless parameters. 

He identified five regions of parameter values vhich were 
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used to exhibit five different types of multiplicities. He 
then applied the same methodology • to the case of variable 
initiator concentration i.e,/ vhen the change in the 
initiator concentration cannot be neglected. However, 

Choi's methodology does not have any provision to include 
the gel effect. In fact, for bulk polymerization systems 
or for solution polymerization systems with fs 0«6 , the 
presence of gel effect cannot be overruled. 

In the first part of the present work, an attempt 
has been made to improve the correlations for the gel effect, 
Ihere are two parameters , ©. and © (See Table II) the 

t p 

13 

model of Chiu et al. which need to be improved in order to 

on batch reactor 

match experimental data In the second part, a study of 
free radical solution polymerization of methyl methacrylate 
in CSTRs has been carried out, with an objective of investi- 
gating the steady state multiplicity patterns of the reactor. 
Ihe influence of both the gel and glass effects has been 

13 

accounted for in the present work, ..using Chiu et al's gel 
effect model with improved © and © correlations obtained 

t p 


in the first part. 



CHAPTER II 


FORI'^ULATION 

line first step in formulation is to identify the 
kinetic mechanism. The kinetics of free radical polymerization 
is summarized in Table I. The reaction mechanism adopted 
here consists of three major steps, characteristic of all 
free radical polymerization processes. They are initiation, 
propagation and termination. In the chain initiation step, 
free radicals are generated by the thermal decomposition of 
an initiator such as AIBN (2, 2' -- azo bisis obutyronitrile) . 
They result in two primary free radicals, R®. In chain 
propagation, primary radicals thus generated quickly react 
with monomer to form long polymer radicals, P^: Propagation 
ends when either a glass is formed (and k^ — ► O) or an 
equilibrium monomer concentration is reached. Termination 
occurs when two radicals react via a bimolecular process. 
Radicals could terminate via two paths. The microradicals 
can combine to form a single dead polymer chain (called 
combination) or form two dead chains (called disproportionation) 
Although radical transfer to solvent and to monomer has to be 
accounted for for accurate prediction of molecular wei^ts, 
in this study, for the sake of sinplicity, these are neglected. 
Also neglected are the branching reactions, thermal initiation 
of monomer, and all other transfer reactions. These are not 
e3q>ected to influence multiplicity characteristics significantly 
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The second step is to develop the species mass balance 
and energy balance equations from the kinetic scheme of 
Table I. A CSTR of volxjme V is considered hece with 
inlet and outlet streams as in Figure 2, The symbols are 
e3<plained in the Nomenclature, Details of the baJLance 
equations are presented in Table II. A little reflection 
on the definitions of the Q and M terms shows that the 
volumetric flow rate of the product stream is given by the 
sura of the contributions of the monomer (ra), polymer (p) 
and solvent (S); 



( 1 ) 

(a) 

(b) 

(c) (2) 


In Eqn, 2 is the gm of monomer/min, being converted 

into polymer, "E = f /(1-f ), and f is the volume fraction 

S3 S 

of solvent in the feed, . It is assumed that the density of 
the solvent is constant, and that the feed consists only of 
monomer, initiator and solvent. 

- The monomer conversion, is defined in Table III. 

Defining the volume contraction parameter, e , as 




C3ts , ,Qp 
P m , l^p > 


Figtare 2 • 


Schematic presentation of the continuous 
flow stirred tank reactor (CSTR) , 
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TABLE II 


SPECIES BALANCE EQUATIONS FOR A CSTR [Fig.2] 


Monomer Balance: 


[ MT ^ “ C' r M T = 
0*^0 ’ 


k 7^ r NT V 
P o 


(a) 


Initiator Balance: 


Q„[ i]^ - Q[i] = K^[il V 


(b) 


Live Radical Balance: 


n = 1 : 


Q [PJ = (2f [I] - k^ [PJ (c) 




n > 2 


“[ 


< ’'b [»J ‘ [Pn-l] - [ ^n] > 


f "2^ ) 

t. L n j o ^ 


V 


(d) 


Energy Balance: 

Q P C_ (T^- T) + (-4 H ) k f M m V - hA (T-T^) 
o mo P f ^ ' j--' p t_ /o . c 

Ce) 


O 
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—p 


(3) 


and with the help of Eqns. (2) and (3) Eqn. (1) can be 
simplified as. 


Q 


Q. 


(1+^) 


(1- S X ) P 
^ 1 mo , 

L 




m 


] (4) 


Tlie various model equations of Table II can be made dimension- 
less with the help of Eqn. (4)» They are presented in Table III. 
Ihe various dimensionless variables and parameters are also 
defined there. 


For the live radical balance, the quasi-steady state 
approximation (QSSA) can be taken for the live polymer. 

Although QSSA for live polymer fails at high conversion for 
nonisothermal batch reactors, several researchers have used 
the same for free radical solution polymerization in CSTRs^^'^^ 
without much error. Ihus equation (d) of Table III reduces 


to 




5 k 


t 


(5) 


Equations (a), (b) of Table III and Eqn. (5) can be combined 
into a single equation. 


F(X^,Da, a, p , 6 ) 


k(T^)2 


9 0 


(1 


Xi) 




k 

r s ^ Da 1 = 0 (6) 

k(Tf) ^ 
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TABLE III 

SPECIE5 BALANCE EQUATIONS IN DIMENSIONLESS FORM 





TABLE III (CONTD. ) 


Where 


= 




O ‘- o 


H - M 
o 


M 


^2 = 


Q o In - Q [IJ 


Q o lI] o 


X, 


T - T, 


a 


h A 


^ k(T.) V 

‘mo p ' 


/8 = 


(- A H^) [ M 
* mo P f 


T — T 
C 


T = 


V 


E = E^+ (Ed-E. )/2 

P ^ 


^(T^) = k° exp (-e;/RT£ 


2£ r I u 

) ( — )% 


(!-• X-) ^ 

S = [ ± 552 +i;]/(l+ii) 



In order to obtain the relation between and in Tabl e III 
equation (a) is multiplied by p and thoi subtracted fram 
equation (c), giving, 

p X- a Da 

X_ = i + (7) 

■ 1+ a Da 1+ a Da 

For isothermal batch r eac tor simulation, mass balance 
equations relevant to the kinetic scheme of Table I, are 
presented in Table IV, using the present set of variables, 

Ihe constitutive equations for the gel and glass 

13 5 

effect models of Chiu et al . and Boss and Laurence are 
shown in Table V, While the Ross and Laurence^ model shows 
a discontinuity in the value of k and k at the onset of 

t p 

the gel and glass effects respectively, the model of Chiu 

13 . 

et al. gives k. and k in the form of continuous functions, 

t p 

and hence these are easy to handle. Termination and propa- 
gation rate constants in the absence of the gel and glass 
effects predictably bear exponential relations to the reactor 

temperature, T. ©. and © also vary exponentially with 

t p 

temperature. Althou^ © does not vary with [ l] , 

13 

indeed a function of [ l] (discussed later). Ihese 
functional relationships are also presented in Table V, 

Eqn. (6) defines a surface in the parameter space 


defined as 
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TABLE IV 

SPECIES BALANCE EQUATIONS FOR ISOTHERMAL BATCH REACTOR 


Monomer Balance; 


dX, 


-5f = k (1-X,) ^ 

dt p 1'^ o 


(a) 


initiator Balance: 


d [ l] 
dt 


[I] 


+ 


■ill 


1 - eX- 


(l-Xi) kp (b) 


Live Radical Balancer 


O = o_ ^ k + 2f k [ I] - k (c) 

dt 1- 6 ^ P ^ ° 
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table V 


CONSTITUTIVE GEL AND GLASS EFFECT EQUATTC^S 


Ross and Laurence Model : 


Gel Effect: 


to 


0.10575 exp [17.15 - .01715 

(T - 273.2) ] for > 0.152 


.-5 


= 0.23x10 exp (75 V^) for £ 


O. 


Glass Effect: 


^ = 1.0 for 0.05 

po 

= 0.71x10”^ exp [ 171.53 V ] 


where = 0.025 -h « (T - T ) + fa - T^) 

f m m L p P 


(T - T ) 1 0 ^ 

m ' m P 


13 

Model of Chiu et al. : 


Gel Effect: 


'to 1+ ®t ^to ^o 


Glass Effect: 


152 



TABLE V (CONTINUED) 


where. 





2.303 0 

c 

r= 

exp 

[ 21 




A+B 0 
m 

9 

P 

= 

9° 

P 

e3(p / ^9P V 

RT ’ 


= 


exp ^ ^©t j 




RT 

k 

Po 

= 

k° 

Po 

exp (-Ip/RT) 


=r 


exp (-Et/RT) 

o 


o 




k° 

exp (-Ed/RT) 



1' 

-^1 


1- 

x^+ e 




X^(l- e ) 


1- 

s x^ +e 



and this surface is called the steady state manifold of the 
system. A bifurcation parameter is said to be the one which 
can give rise to transitions from unique to multiple solutions 
or vice versa in a plot of conversion vs. this parameter / 
for a given set of values of the other parameters. For 
example/ dimensionless mean residence time. Da is usded 
as a bifurcation parameter in Figure 1. Plot (e) of Figure 1 
e:<hibits four bifurcation points (transition point between 
uniqueness and multiplicity) in the transition scheme of 

19 

solutions. According to Golubitsky and Keyfitz , 

★ 

the surface in the parameter space, p at which the continuous 
change of a parameter causes the appearance or disappearance 
of an S or inverse S pattern (Figure 1(a) and (b) respectively) 
is called the "hysteresis" variety. The following conditions 
are to be satisfied by the hysteresis variety^*^: 

(a) 

(b) 

(G) (9) 

Similarly, the isola variety defines a surface in the parameter 

i( 

space, p , where the continous change of one of the parameters 
causes the apprearahce or disappearance of an isolated branch 
in the conversion vs, mean residence time plot, and satisfies 


F ( Xw Da, p ) = o 

± SX 


5 F • 

r- — ( X., Da, p ) = O 


2 

( X., Da, p*) = o 

' 5^2 1 V, 



2 2 


Eqns. (9a)/ (9b) and 

Us Da, p ) = o (10) 

Introduction of the influence of the gel effect on the 
propagation and termination rate constants makes the mass 
and energy balance equations very complex and consequently 
eqn. 6. For this reason, it is extremely difficult to 
eliminate Da from eqns. (9a) - (9c) to get an explicit 
relation between a and p for the generation of an a - jS 
plot. Function F, in eqn. (6) is differentiated twice with 
respect to to give three non-linear algebraic equations 
in X^, a i’p / Da and 5 , which can be solved simultaneously 
with the NAG subroutine C05NBF to give values of X^, q. and p 
Da for given values of p- and 6 , This gives the a-p plot 
for the hysterisis variety. Similarly, for the isola variety, 
function F is differentiated once with respect to Da and the 
resulting equation is solved simultaneously with eqns. (9a) 
and (9b) to give an a - p plot for the isola variety for a 
given value of 6 , ihe whole oc- p plot thus generated is 
shown to be divided into various regions, each representing 
different multiplicity patterns of the conversion versus time 
plot. Particularly, the whole a- p plot has been divided 
into five major domains. Each type of multiplicity pattern 
can be verified by picking up parameter values from tliat 
region and generating the conversion vs. Da plot. For batch 
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reactor simulation and for the study of steady state 
multiplicity of cSTRs, Pmma polymerization in CSTRs is 
taken as an example and the various physical properties 
and parameter values for this system are presented in 
Table VI. 



TABLE VI 


VALUES OF PARAMETERS FOR PMMA 


14/ 31/ 32 


£ 

= 

0.58 (AIBN) 

A 

= 

0.168 - 8.17x10"^ (T-387)^ (Ref, 

B 

= 

0.03 (Ref. 13) 

s 

= 

0,4 cal/g °C 


= 

0.973 - 0.001164 (T-273) g/cm^ 


= • 

1.2 g/cm^ (Ref. 13) 


= 

30.66 kcal/mol (Ref. 13) 


= 

4,35 kcal/tnol (Ref. 13 > 


= 

0.701 kcal/mol (Ref. 13) 

^©p 

= 

27.82 kcal/mol 

^©t 

= 

34. 176 kcal/mol 


1 = 

13.5 kcal/mol (Ref. 35) 

k° 

^d 

= 

6.32x10^^ min"^ (Ref. 13) 


= 

7 

2.95x10 lit/mol-min 

Pq 



ko 

to 

= 

5.38x10^ lit/mol-min 

e 


— 16 

5.48x10 min 

p 



< 

= 

— 22 

3,372x10 mol-min/lit 

0 

z= 

r ■ 

™ - 0.1946 + 0.916X10- ( 

iv 

°^p 

= 

0.00048 (Ref. 5) 


= 

0.001 (Ref. 5) 

T 

P 

= 

387.0 (Ref. 5) 

T 

m 

= 

167,0 (Ref. 5) 



CHAPTER III 


RESULTS AI<ID DISCUSSION 


It has been mentioned earlier/ that we have used the 

13 5 

gel effect models of Chiu et al. and Ross and Laurence , 

13 

in our work. Results using the correlation of Chiu et al . 

are presented first* A carefxil survey of the various 

3l 3 36 *“ 39 4 1 

publications of soong and coworkers * * reveals that 

o o 

they have used different values for the parameters ©^, ©p/ 

^© ^© in their correlations for ©. and © (Table V), 

t p t p 

13 . . 

Chiu et al, had matched experimented, conversion - time 

7 

data on batch reactors of Marten and Hamielec at three 

temperatures, 50®C, 70"C and 90* C quite well, using values 

of © and ©^ shown in Table Vll. Ihe agreement at one 
p t 

temperature/ namely at 90*0 for two different values of [ Ijq 

13 

is shown in Figure 3 (curves a) Chiu et al. , tried to 
correlate the three sets of values of © and ©. given in 

p -t 

Table VII with the exponential relationships; 

(a) 

(b> 


© = ©° exp ( Eq /RT) 

. ^ P 

©t = ®t ^ ^© 


( 11 ) 
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TABLE VII 


9. and 9 Values Used 
t p 

Experimental Data of 

15 

by Chiu et al. 

7 

M ar ten -H ami elec . 

in Fitting 


Initiator loading# 

[ I ] o 

Polymerization 

Temperature# t®C 

^t min 

A 

9 

p# min 

0.0258 ' 

50 

1500.00 

3500.00 


70 

49.00 

250.00 


90 

3.80 

3o4oo 

0.01548 

50 

2330.00 

3500.00 

• 

70 

83.00 

250.00 


90 

6.30 

30.00 



0 5 10 15 20 


t (min) 


Fipjre 3 ; Comparison of batch reactor simulation 

■ ri'i-'nil hs fur PllilA ix">l yinf7"i. tJ.cjn u.'i.i.n n t ''’J. jmi' J I'-f 

^ r and E from Table VIII 

t p 

(column £>) -Vi th experimental data(o,« ) at SC^C 
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IhuS/ plots of In (1/9) vs, 1/T would be strai^t lines, 
which would give the values of the parameters 9^, 9^, 

and E Iheir values are reported in Table VIII in 

fey , fey • 

column (b) . When these values of the parameters are 

used at 50°C, 10°C and 90°C, it is observed that the 

calculated values of 9^ and 9 are quite different from 

t P 

those given in Table VII. a comparison at 90®C is also 

shown in Table VIII. Obviously, conversion “ time plots 

do not match experimental isothermal batch reactor data 

at different temperatures, as shown in Figure 3 (cuirves b) 

37 

for 90°C, Louie at al, used values of the four parameters 
0°, 9°, , and E^ as given in column d. Table VIII. 

t p 

using these values it is observed from Figure 4 (curves d) 

that the conversion vs. time plots at constant temperature 

again do not match experimental results too well. Another set of 

40 

parameter values were obtained by Kapur at al . (Column c. 

Table VIII) to match nonisothermal plug^-f low reactor simulation 

36 38 

results of Baillagou and Soong ' . The corresponding 

results for isothermal batch reactor polymerization at 90 *C 
are shown as curve c in Figure 4, Poor agreement with 
experimental results is once again observed. In the present 
work, we have used the values shown in column e. Table Vlii, 

since these match constant temperature simulation results 

41 . 

of Carrat et al. on CSTRs. ihe values shown in column b. 

Table VIII fail to duplicate the results in Ref. 41' . 
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TABLE VIII 


Comparison of the Values of the Parameters and 

^ ®t p 


Used by Various Workers 


13,36-39,41 


Curves a 

Parameters 




r II = 0.0258 
'■ mol/l) 


ri] ^=0.01548 
° mol/l ) 


G (90“C) 
P 


Ref No, 



5.48x10“^^ 

5. 48 14x10 "■ 

16 

5.48 14xl0"" 

■16 

5.48x10“^^ 

-> 

-22 

1.14x10 

1. 155 3x10'’ 

'22 

1.1353xl0" 

.22 

— 22 

3.372x10 

- 

17400 

17200 


17420 


17200 

- 

14000 

114100 


13982 


14000 

80 

2.9021 

1.695 


3.05 38 


4.9478 

30 

4.836 

2.8253 


5.09 


8 . 246 2 

hOO 

30,79142 

40.567 


29.31 


30.79 14 

13 

41 

40 


37 

Present work 
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Figiire 3 (curves e) sho\«that this choice of parameters 
again leads to poor agreement with isothermal batch 
reactor data# 

CSTR simulation results have been generated using 

the parameter values of column 6/ Table VIII, at constant temperatures 

of 50®C, 60®C, 70®C, 80®C and 90®C. Ihese are shown in 

Figure 5. Also shown in this figure are the simulation 

41 

results of Carrat et al, using, the parameter values in 
column b. Table VIII. There is some difference observed 
between the two sets of results, particularly at low 
temperatures, but this was the best agreement we could 
obtain without using an elaborate optimization algorithm 
for obtaining the parameters. These parameters (Column e. 

Table VIII) have been used finally to study multiplicity 
of steady state in CSTRs. 

Figure 6 shows some constant temperature simulation 
results on CSTRs using the parameters in columns a and e 
in Table Vlll. Considerable disagreement exists, particularly 
at the lower temperatures between these results. It is 
known that 9^ is a function of the molecular wei^t. This, 
in turn, is a function of the initiator concentration in the 
reactor. Chiu et al . , used 0^ as a function of C 1 ] ^ the 

initiator concentration in the feed (to the batch reactor) . 

However, we replaced this by © = [11 ), and accordingly. 
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generated simulation resiolts on C5TRs. ihese are shown 
in Figure 7, We observed that this does not lead to very 
different results compared to when [ I used^ and 

so to use the latter for multiplicity studies. 

In summary, we observe that there is quite a 

confusion regarding the parameter values for and ©^ 

in the literature, and even the same group of workers have 

used different parameter values in different publications. 

Our own set e (Table Vlil) was thus used for multiplicity 

' studies since this gave simulation results on constant 

41 

temperature CSTRs, matching those of Carrat et al . , 

13 

Along with Chiu et al. *s model, the model of Ross 

and Laurence^ for the gel and glass effects is studied 

extensively to predict its applicability for free radical 

polymerization in batch reactors 

and CSTRs. For PMMA polymerization, 

values of the various parameters in the model of Boss and 

Laurence are shown in Table IV and the model is presented 

5 

in Table V. The Model of Ross and Laurence for k,/k^ and 

t t,o 

To 

k /k is used along with the k. , k values of Chiu et ali 

p po to po 

to generate batch reactor simulation results. Ihese are 
shown in Figure 9 . Results for varying initiator efficiency, 
f, as Well as a constant f of 0.58 are shown. Ihe results 
using variable f are close to the experimental data. 



T(min ) 


Figure 5 . 


c 

CSTR simulation results with parameters 0^ 
0°, and from Table VIII (colum'. e) 





0 ( 


I 


Figure 6 


T (min) 

• CSTR simulation results using parameters, 
0°/ E_ and E_ £i?om Table VIII (column a 



GSTR simxiLation results are generated for constant £, using the 

Ross and Laurence equations for the gel and glass 

effects. these are plotted in Figure 10, Ihe Ross and 

Laurenace^ plots show two discontinuities in each of the 

curves, one at the onset of the gel effect and the other 

at the onset of the glass effect. Ihese are associated 

with the use of different equations for the gel and glass 

effects in Table V before and after a critical value of the 

free volume. Ihe results using the parameter values of 

Table VIII (column e) at 50*C, 70®C and 90®C are also shown 

in Figure 10 for comparison. Althou^ they match well at 

90“C, for other temperatures they do not match. Obviously 

multiplicity features are ej^ected to differ significantly 

for the two models. In this work, we have used the Chiu 
13 

et cQ. , model for the gel and glass effects with parameter 
values as in column e of Table Viii to study multiplicity 
of steady states in CSTRs. 

35 

Choi studied the steady state multiplicity patterns 
of CSTRs without the gel effect, considering free radical 
solution polymerization of PMMA with fs > 0,6, and generated 
an ct-P. plot to depict the regions of various multiplicity 

20 

patterns following the technique of Balakotaiah and Luss 
However, he did not consider, in his methodology, the change 
in initiator concentration due to reaction, i,e., [ I] = [ i]q* 
In the present study, the change in the initiator concentration 


Figure 7 


C5TR simulation results with pa^aTieters 0 

^0 Table VIII (column e) 

■t 'P 


for © = © ( fl] ) 




TrCmin ) 


Figxire 8 


Comparison of C5TR simulation results with 
o o 

parameters 9 , 0 , E and from 

t p 

Table VIII (column e) with the results 

41 

reporduced from Carrat et al.'s plot. 




[l]o=0.02W I / 

/ /yy 

y / 

y y y 

y / ^ 

yyc 

yy 


f 

/ 0 . 


015 48 


[Ilog.mol /lit 




t (mm) 


Figure 9 : Batch reactor simulation results using Boss 

and Laurence gel effect model. Results for 
varying, f =0.487-2,51 (dotted line) and 


£ = 0,58 are presented (solid line) 
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has been taken into account. Figure 11 shows the ct-p plot 

for the hysteresis and isola varieties. In this figure/ 

35 

curve a shows the results of Choi for the hysteresis 
variety in the absence of the gel effect. The rate 
constants are different than those given in Table vi. 
and are: 

k, = 3.57x10 e min 

d- 

k = 4.2x10® e “^300/RT iit/(mol-min) 
po 

k. 2.112X10^^ g-2800/RT iit/(mol-min ) (11) 

Curve b shows the hysteresis variety using the rate constants 

13 

(in absence of the gel effect) of Chiu et al. / as given 
in Table VI. in addition, the assumption [l]= tlJp 
not used, nor are density changes accounted for. similarly 
curves c and d give the corresponding isola variety plots 
for the same situations. 

Simulation results for CSTRs for some values of a 
and p in different regions of Figure 11 are shovm (Fig, 12) using 
the r ate constants (in absence of gel effect) of Table VI 
and not using [l ] = [I ]q* Uniqueness, hysteresis as 

well as isolas are observed under appropriate conditions. 

The last part of the present work is devoted to 
the study of steady state multiplity with the gel and 
glass effects for PMMA polymerization. The functions F, 




Ftoss and Laurence gel effect model with 


Chiu et al.'s model. Parameters ©°, 

t p ' 

and E_ are taken from Table Vlll 

t p 


(colxjmn e), f = 0.58 


41 


^ and F_ are so stiff that simultaneous solution 
of the two sets of eqns, (9a)--(9c) and (9a), (9b) and 10 
is very difficult. Ihese simultaneous equations are very 
sensitive to changes in the initial guesses using the NAG 
subroutine C05NBF. For this reason and for lack of time, 
an attempt has been made to obtain an approximate a - p 
plot from simulation results, rather than from the technique 
of Balakotaiah and Luss^^ 'Hiase results are shown in Figure 13. 
Simulation results taicLng some of the parameter values from 
Figure 13 are presented in Figure 14. These plots have been 
generated using data in Table VI and column e of Table VIII, 
It is hoped that precise boundaries in the plot in 

the presence of the gel effect will be obtained in the near 
future. 
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0 10 20 30 40 50 

rChrs) 


Figure 12 ; 


CSTR simulation results for free radical 
solution polymerization of PMMA using 
different values of the parameters 


and P 
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igiire 14 .-simulation results with gel effect for 
solution polymerization of pmj/a in C5TRs 
using different values of a 


and g 



CONCLUSION 


Modified parameter values for the gel effect 
equations of Chiu et al . give closest match to the 
CSTR simulation results of Carrat et al. and can be 
used for the stiidy of steady state multirlicity of 
CSTRs with gel effect. Although only two types of steady 
state patterns are obtained from the simulation results for 
the case with gel effect a specific bifurcation diagram 
is needed for the possibility of other types of steady 
state patterns in this case. A knowledge of stability 
of these multiple solutions may help efficient control 


of the reactor. 
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pmk PotmM zhm>i t'J tbir 
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4 f 8 s 

?^LPH/i = 

BETA = 

OZhfk = 

« m 

5/00 = 

EFF s 


POOJifA, •sDOLLrrbE' ::y«REOrtriD,si 


EFFSC.t i’O Ttlftr 


/ R 


l^ARtABLCS 

parameter df Tiis: 

FDR CO'/OO), 

REA.r transfer parameter 

TEAT 3S,N£.RA.ri3^i' PARA'iSIER 
0IM2|ilSi.3M(^ESS ZOOIiA.MJ" TEMPERATUR 
0/03 

RArrD OF orFFJ-STYriiss ^?ith gel 
<^ rr,-i G£L' EFFS'^r.t 
IHiriArOR EfErfHlnOT 

ACXr.VAflOM Ef^S’R'3.i FOR: PROPAGAriO.M' / R 
AGTrvAriOM SMS:RG.jf for- I.fIiriATDR oe^O^P 
AGflVATlO^ EriS:R'G.f FOR' TERMI^ATiO^i / R 
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ABSENCE: OF GELt SFFSCT'^ FOR XMOEf^l = 2 
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PRfi, - SiXPOMENTXALl ?A:RT OF K.:PO 
LI if SI RADiGALi ?D:?.JLA;ri,OM 

pRootJCi OF loo; a:no' kp. 
temperatlirs. ,n oss..' X:,* 
fee:o' rs;MpERAXJ.Rs; 

OAMXDUbGR MUM3E:R1 
vjQMOMSR (/OLtJMS: FRA;CirX.OM 

DEMSITX “■ 

^DESSlXi 
ifOLJME. 

ROMO / ROM 

DEMSITX RATIO R'OM/ROMO 
MOMOME.R CDMi/ERSTDM 


ET ■ : 

^ - 

ISDEX2 

13 ^ 

IMOEXl 

KO 

K.r 
K r 3 
KTOO 
K? 

KP3 
KPOO 
L03 
LOOK? 

T 

TF 
DA 

PRIM 
ROM 

■■ 

SPSS 
F3- 

DR' ■ 

XI 
X3; 

m Kp :m m m m tm m'm m 

DTMENSIDN' X(2D3 
INTEGER 1 T M A X , r £ R # M f N S I G , I M OSX I IN oex 2 > I « T 
COMMON alpha, BETA, TA. a.,I,M5exi'fi:MOeiX21,FS, rf 

external ryM 
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iFFECX' 
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F WjlB* < 






FS = 0.6 
TF = 300,0 
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DSL«H’ 


13. 5 
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7^9 


5 56 ' 


581 

533 


1531 

1502 


18} 

555 


6 0 

131 


5DC0.0 

RO'IO = 0.973 - 0.001164 ♦ C f f *i 2:7 3.. 03 

S'lfeTF = (1,0 -FS3 * R3M0 ♦ 1 0 30.,0/lSMS;:4' 

beta = Def,HR ^ E 3 ,SF/(R 0 CP * Vf) t TifPg- =I',B.ETA. 

i»fRl,TE(22,70a3 ALP8A,3£fA 

3')Cr 'aEirAs.'rB-U.'O/)^, 

rAu=o 

IFdi^OEjfl.E 
I,P(IN 0 EX 1 .E__ 

DO 102 M=l,ITr 
rypE 556 
F0R«Ar€5X,'r 
XCi)~».957 
XC 23 = »/.B 5 
X( 3) =0,613 
X(4)=0.291 
X(5)=0.1 
K(6) = 0,05 
IFdODExi.ED 
IFCI^DEXI .EP 
rAUrTAU-dlMC 
IFCIiiDCXl.ED.'l) 

' "d'l) 


ii) I'irsloo 

2 3 X?IT =. 15 0 


hi 3 SR£.*> 


'13 dNC =3 0 
'2 3 I’r.NE = 1 


IFCIMDEXl.E 


dME=.TAiJ/. 5 i 3 
If PS 5 8 0 ,ri^'Si 


IFd^'DtXl,E^.23 riME=,TStU 

i:f(T6oex1,e 5.'23 rxps 531 , rim? 

FORMAT! 5X, dlMS l.'l- '4l,^^U^£S =- '>ST-3','6) 



CAii£"iREALrcFrfM7i,Ps7iPS2"irI:',iiiE7i7x7i5iAXd'R'8) 
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DO 555 J=J,5 ^ 

IF(X{J3,GT.1,03 3.0TO 555 
IF(X(0),LT,0,03 S'OTO 555 
Xl •= X( J3 
rXPE 181, XI 

WRITE (22, 181) Xi^ . . , . 

FORMAKBX, 'SOLUTION - *^E13.50 

C.ON'fifjUE 

WRITE(22) 

comtinue 

F3RMAT(5X, ' DfPE: I.N' T8E ^ALCl.g: Of. INOSKl V) 
FORXATISX, ' d'PE., l.N- TNE. 7ALUl! Of, 1:105X2 V j 
STOP 

END ' 

function FUNCxI) ' 
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rifPE^=rKF 
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rm fw'h f rr M. m ■w jLf ri 'fltt m ft it '*»”»» T » ' 


EVftLUftriwIN OF the: FUMEllO'^ , f 

m m'mmm.m mm mm «»W m mjm m'm'"* "* 

J.F’cr>M0FXl,E0.2)'S3rb I'u 
KfsKTO ! KP=XP0 

3 r I rt ij F 

?J7CJ=2r0=»*EFFl'1.01‘CKPl'<'2)l=K0/IK:jrr'CfW*23 
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OPlilM = ((1,3 -; XI) * JEeSK -! l..<3 f SPSS ♦ Kl -i BSi/ 
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n = -THETA?. * X.P0 * 3C3Xi-/i l=*2* + CrHETAP * DXP002 + xpo A. 
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SPSS = 0,1943883 ¥ D .-9 1 & 5 SE-B* ( T -: 273,0 ) 
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I3/2X , ' U 1 5 .1. ,Kt 3 i , 3Ei 3 . 6 ) 
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■P\FBn> * V i 4 \ /■/•Tfit Jx»'iT ”1 r» if - 



£,\[ H.S v/ftLiyS 3; 


3.5S S s:S^ / 5 X . ' X C 1 3 » X £ 2 j r a ( i ) =.' / j 
1,^ 3'SX V 5 X *. * I M p Ci =4 [J\ S ^ C rt 3 i. 


15 RAT^ ^SmsTaItS V5X, 'I^DEX=2' '«:2A;?JS SOOfISS RME' 
IFCBSTA.GT.l.O) SfPP 
33-rO 8 31 
SP'PP 
S*ID 


3jdKjijj'r*yE Fm hf; E^f6b.aArj3M- trs FUMCjiaMS' 




SiJtiRailfT'iF TZ i (7, K, F^SC, IJu'AGp 
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INTEGER 'iZU, , T-43i.X , f Fu'Aw., S , , 
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DEFIMrin’ff OF PA?4,<nHL€R'‘ MUH'SSR ^ , , . . / ' ..... 

« »;■,« » «U • m'm m'm m^m 

E = EP + Q5f (E3i- Bfl 

KF =’KP00 * EXPC-E/TF) ♦ Sp'RSTl- 2S.0 * E.FP. ♦ 10 * KOO/MOO J 
TAO = DA/.KF 
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DEFI^'irrON 
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EP 
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FT 
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KPOD 
LDO 
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TF 

DA 
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EPSS 
F3 
DR 
T At! 
fHETAP 

XI 

Kl 

DCOXl 
DADXi 
D2 ADXl 
DPC,DXl 
DfDXl 
DBF 
D31’2 
> 

OLDFP 

DLDD 2 

DLOKP 2 

oxroxi 

DXrUDX 

PFPDX.1 

DKPDDX 

DKOOKi 


OF VARIABLES 

- PARAMSrSR^ JF THg: FOOJITA DOLLFiLsi; ZUrREL a.T LD i 

FOR 0 / 00 ). 

= HEAT transfer PARAl'^ST'ER 
= HEAT 3 S;i 2 HATi 3 )J- ?,ARA«!SrER 
= DIHSNSLOMLSSS :D 3 L!A.Yf TEMPERArOR i; 

= 0/0 3 

= RAT to OF 0 IFrJSl;/TTt.-S ,5 «l fH GEL* EFFECT' 10 TiiAl 
/^ir-j: G.ELI SFFOF'.t 

= t.^I.riAfOR EFFtClISHCif 

= ACTI.Va^l 0 ^i FOR: pROPAGATlO.i /■ R 

= ACi't.YAnOrJ EMF'RG,1 FOR INI FI ATOR lECOap. 
s ACriVAriOSJ ESiFRG.f FOR rSRHlVAnOv / R 
= SOLOENr YObOMS: FRACiriO^ l«i THE FEED' 

= PARAs^ETER INO t CiATT G. C'HOI ' S RATE ClNsiU..^! For 

INDEX =• I A^D SOOYG.'S RATE COMSIA'^T FOR IMD 2 .X = 2 
= INi.l’IArOR COMCSlNtRAliOM r GM - vsjf,/Li.r 
= PAKAMETER I SO ICiATirN G- GEL EFFECT FOR ^GiiU = 1 a^o 
arsemce: of GEL! effect for NG-EL =2 
= INiriATOR OECDMPiOSlTlON RATE CO^SFA^tf 

= PRE. - £;xpo\ieNri;A: 4 i part' of xd 
= rsRTiHArro'J RA:rs: constant 

= ^ERv!I,^f^lTiOM rate; COVSTANT IM ABS'SHCE of >»ey.' tFFECiT 
S ?,RS - E.XPO.'^SNPI.ALI PART OF KTD 
S »R:)PaGA,T10Y RATS: const »A. iU’ 

= propagation rats; CrONSTANT IH absence of glass EFFEC.r 
= Pftd. - EXPONENTTALI PART OF KPO 
s LT Y 51 RADICAL ?O;?.JLA.Ti:0i'i 
= PRODUCT OF LDO. AiNO' KP 
= TEAPERATJ.RS. IN OSG..f Xv 
= FESO' TEMPERATJRS: 

= OA'^XOHLlER NO'^SSRi 
s MONOMER' V0WIH5! FRA:C]r’I'ON 

= DSNSTTY OF THS: HONOH-S'.R AT THE REACTOR TS.'^ PER AT ORE, 

= OENSITY OF THS> OTNOMSR AT THE. FEED TEMPERATURE 
= NOLJ.MS, CONTRACiriON PARAMETER , ( ROp - RJMJ/HOP. 
s RO^O / ROM 

S OEN'SIT/ ratio , R 0 :M:/RD«.D 
a mean. R 5 SI, 0 S.NCSi TtME! I,N Mp.* 

THETAT =. PARAMETER'S LN. THE. GLASS AND' G'FL EFFECT' MOOfiLi 
OF CHIU ST*. < Alii. I R'ESPECTlVEuiY 



IKC. 

FTKST OeRINATIXS: OF, Ci NIT'H RESP-SCT TO Al 
FT'BST OeRIYATl.vg! Of', a: to XI . 

s SE'w.OiND OSRLVAr'LN'S! OF A «,R,. TO. XL 

a Second osRLYArtytl of. c: «.r.' td xl 

a FfRST OERINATTNS; OF TEM,PS,RATURE i' N.R.t TO Xi 
a first DSRI.NAriYSl OF U/T) W.R.' TO ki 

a SECOND DERI-VAtiySl OF U/D ?f.R-! TO X» 

a first OERINATITS; of loo TO Xi 

a FIRST O'ERlYArrYS; OF LOOX§ W.R'.t TO kI 
a SECOND DSRI.VAniNSl OF LDD »i,R.- TO XL 

a SECOND OSRr.NArilYSI OF LOOK? . «.a,'T-i Xl 
a FIRST OCRINAII'N’E: OF ...XT N,B.. TO XI 

a FIRST OERINATrVS; OF KTO ^.R.- TJ Kl 
a first OERINATLsIS: of XiP. N,.R,' TO XI 
a FIRS'T OSRI'NATIlfS; OF K.P .0 W.R.- TJ XL 
a first DERINATirS' of XO- T'O XL 
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BUCH REATCJR S£Ma:,ftl’£,a»| ; RXVi:3,£= «1 KJfTR 


Dl .'^ ENSlQt ^ iPRI ' lEClO ) » 
IMI’KGEH riimgs. 

REAIj KP , K£ ,<D , 

i RPKP 3 , IfKrJ 

3 PEN tu.ur = 22 , OZ'.sllZi: 


a IB ) 

SP 3 , mi , £3 
'i 35 i <: ' 


, , LJO f A 

■ # FTuE = ' 


H , = P.OOi 
fR'Ep = 1000 
K = 0,0 

)Cs:ni.> = 33.0 

lEOlJwX = 1 

= 33D.0 
=■ 363.0 


TF 


13 


u, 01548 


1.4 lTI-%0 cn -'4011134 3 


yen = a,u 

{(2) ~ I a 
7(33 = 0.0 


111 


4 Rl . TE ( 22 r 100 ) 
44R1 1’EC/?., 203) 
. 4 RirE ( 22 , 303 ) 

K = R-ONGE £3 . 
IF ( K , 4 F , U 
XI s ni ) ; I 


10 


7 PHi : s ^£ 

ij \J 113 


= 7(2) 


# X r 1 ) 
bOO S. 1(3) 


XfO = 5,88 E9 ♦ EXPC - 3 52,7 9-3«/ir ) 
KP9 = 2.95 E; * 2XP( -2190.71^:1') 


KO 


6.32 Et6 * 


XPC -:i5 430-.0-7ir; ) 

1- 0.016567 2;-13. f ( T - 273.0) 
17 S.-06 * (f':-) 3a;,0) I'* 2 


EPSS = 0.1945833 
EPSS = - EPSS 

k = a. 168 - 8 

THETAP £-16 ♦ EXPt IfSB 2 .0/;! ) ^ ' 

r'iSTAX = 1.1 353 S-22/.I0 * ex?,C 17120.0 /!) 


9Hin e (1.0 -I Kl)/:(1,0 + EPS'S'^l'^Xt ) , ^ 
: = EXP( 2.303 * PHiX/tA + 3; PRIX)) 


XpKPt ? = 1 , 0 / { 
KrXTO S 1 , 0 /t 


.0 

.0 


KP 

XT ' 


KPKPO * 
KXKro * 


1 

EPO 

KTO 


THSiTAP 

rBSTft.X 


K ’ PiO ' 

xiro : 


LOO / E ) 

( iOO / E ) 


EFF = 0,58 
ypRI*^E (13 
TpRI^iE ( 2 ) 
1 


= KP ' ♦ Cl.O 
= -ED * 1 

.(i.O + EPSS 


- ¥• liO # ^ 
1.0 -1 XI ) ♦ 
XI) 


bOD * EPSS * t/ 


oo-r ') 



pRIhE (33 = f SfF * KD ^ l; -! U * CLDu 2> - * Cl. -3 

1 “ ^J3 * E?S3 ♦! CiiQO ♦t2)/Cl,D ♦ e?S3 * ki: j 

Giro Ul 

il3 JFCX. uE, XEMDJ 3i3r3 116 

srap 

lU I^DONT = .TCDjMf' + 1 

IFCICuilNr . L.<r.f C.REQ3 SOTD l-M. 
rSfPE t f ^ ^ I » 603 

4Rir)if22 , 433) X , XL 

St 

3310 113 


133 

233 


333 

433 


F3P=''U\T( lOX, /BAr:H REACr3R 51 'IJuCr I 3^1 V/3 

= ;fn.6' 3E.G. K.,^ V.iaX.'lMtl’iftlDR 
1 UOADIHG =: 'Flu. 5 ' «3[.y|i. * /yr) 


F3R;'iAr{.i!}x, ' rT«s: i:'j mi’I.' 
f 3R '3 A 1' C 20X # F 1 0 .6 , 1 0 X , £ 1 3 . 6 ) 


, 5X, 'MOND^ER C'3Si/ , 
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1 

■ 2 ' 

2:2' 


,3i 

33'^ 

I 

44 


5; 

66 


SMO 

syaHDurTi'ig RjM'se -; Kurra 

F3,SCTI0k RU«r,£ C 4^ , X , F - X , -1) 

IwrEGER RUM3S r # <v f -«./ 

DI:^ENS1.0-m PHfCSO) , SA2EY(533 ,, srcx) , rCM) 
DATA «/0/ 

M = M f 1 

G3r0 ( 1 , 2 , 3 , 4 r 5) , 

RyNGE = i 
RSl’UKN 

D3 22 J = 1 , H 
S A S/EX C J ) =. X (JI) 

PHi(J) s FtJ) 

XCJ.) = SAVEX(J) t 0,5 4 « * FCJi) 

K = X + 3,5 4 tl. 

R3PIGE = 1 
RETURN 

03 33 J = 1, N 

PHK J ) = NlC J, ) f 2,0 ♦ f ( Ji) 

X(J3 = SAVSX(J) 4 0,5 * « ? ( Ji) 

RUNGE =1 
RETURN 

1)3 44 J =. l.sj 

PRItJ) = PHI(J) + 2.0 * r (0?) 

XCJ) = SAVEXC J.) 4 H >1^ FC J) 

X =. X 4 0,5 4 A. 

RU.MGE = 1 

DO 55 0 a I , N 

XCJ.) = SAVEXCJ) 4 C PHI.CJi) 4 FCJ’)) 4 H26 
J s. 0 



R[JMGE = 0 
RETIJRW 



